function U = poisson_solver_2d(mesh, pde)
A1 = assemble_matrix_2d(mesh, pde, "dx", "dx");
A2 = assemble_matrix_2d(mesh, pde, "dy", "dy");
A = A1 + A2;
b = assemble_vector_2d(mesh, pde);
[A, b] = process_boundary_2d(A, b, mesh, pde);
% full(A)
% full(b)
U = A\b;
end

%% assemble_matrix_2d
function matrix_A = assemble_matrix_2d(mesh, pde, trial_type, test_type)
dNdx = reference_basis_2d(mesh.basis_type, "dx");
dNdy = reference_basis_2d(mesh.basis_type, "dy");
phi = cell(mesh.Nlb,1);
psi = cell(mesh.Nlb,1);
matrix_A = sparse(mesh.Nb, mesh.Nb);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
    switch trial_type
        case "dx"
            for i = 1:mesh.Nlb
                phi{i} = @(x,y) ((V(3,2)-V(1,2))*dNdx{i}(x_loc(x,y),y_loc(x,y)) - (V(2,2)-V(1,2))*dNdy{i}(x_loc(x,y),y_loc(x,y)))/J_det;
            end
        case "dy"
            for i = 1:mesh.Nlb
                phi{i} = @(x,y) ((V(2,1)-V(1,1))*dNdy{i}(x_loc(x,y),y_loc(x,y)) - (V(3,1)-V(1,1))*dNdx{i}(x_loc(x,y),y_loc(x,y)))/J_det;
            end
    end
    switch test_type
        case "dx"
            for i = 1:mesh.Nlb
                psi{i} = @(x,y) ((V(3,2)-V(1,2))*dNdx{i}(x_loc(x,y),y_loc(x,y)) - (V(2,2)-V(1,2))*dNdy{i}(x_loc(x,y),y_loc(x,y)))/J_det;
            end
        case "dy"
            for i = 1:mesh.Nlb
                psi{i} = @(x,y) ((V(2,1)-V(1,1))*dNdy{i}(x_loc(x,y),y_loc(x,y)) - (V(3,1)-V(1,1))*dNdx{i}(x_loc(x,y),y_loc(x,y)))/J_det;
            end
    end
    for i = 1:mesh.Nlb
        for j = 1:mesh.Nlb
            f = @(x,y) pde.cor(x,y).*phi{i}(x,y).*psi{j}(x,y);
            r = gauss_int2d_tri(f, V);
            matrix_A(mesh.elem(n,i),mesh.elem(n,j)) = matrix_A(mesh.elem(n,i),mesh.elem(n,j)) + r;
        end
    end
end
end

%% assemble_vector_2d
function vector_b = assemble_vector_2d(mesh, pde)
vector_b = zeros(mesh.Nb,1);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    J_det = abs((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
    psi = reference_basis_2d(mesh.basis_type);
    for j = 1:mesh.Nlb
        f = @(x,y) pde.f(x,y).*psi{j}(x_loc(x,y),y_loc(x,y));
        r = gauss_int2d_tri(f, V);
        vector_b(mesh.elem(n,j),1) = vector_b(mesh.elem(n,j),1) + r;
    end
end
end